---
title: "Interval Truth Model"
subtitle: "Visualizations of Priors in the ITM"
author:
- name: Matthias Kloft
orcid: 0000-0003-1845-6957
affiliations: University of Marburg
- name: Björn S. Siepe
orcid: 0000-0002-9558-4648
affiliations: University of Marburg
- name: Daniel W. Heck
orcid: 0000-0002-6302-9252
affiliations: University of Marburg
date: "`r Sys.Date()`"
format:
html:
toc: true
number-sections: true
theme: cosmo
code-fold: true
code-tools: true
code-summary: "Show the code"
fig-width: 7
fig-height: 4.5
embed-resources: true
execute:
message: false
warning: false
---
```{r setup, include = FALSE}
# Libraries
packages <- c(
"tidyverse",
"yardstick",
"readxl",
"cmdstanr",
"bayesplot",
"posterior",
"bayestestR",
"rmarkdown",
"cowplot",
"svglite",
"psych",
"here",
"ggpubr"
)
if (!require("pacman")) install.packages("pacman")
pacman::p_load(packages, update = F, character.only = T)
# default chunk options
knitr::opts_chunk$set(
fig.height = 7,
fig.width = 10,
include = TRUE,
message = FALSE,
warning = FALSE
)
source(here("src", "00_functions.R"))
set.seed(35032)
```
# Latent Appraisal
## True Intervals: Informative Beta Prior on Marginal Locations and Widths
### Sample From Priors
```{r}
# priors on marginal locations and widths
df <- data.frame (Tr_wid_splx = rbeta (1e5 , 1.2 , 3 ))
# multiplicative shifting parameters
s <- rbeta (1e5 ,1 ,1 )
# lower bounds
df$ Tr_L_splx <-
(1 - df$ Tr_wid_splx) * s
# transform from simplex to bivariate normal
Tr_bvn <- simplex_to_bvn (
cbind (
df$ Tr_L_splx,
df$ Tr_wid_splx,
df$ Tr_L_splx + df$ Tr_wid_splx))
df <- cbind (df, Tr_loc = Tr_bvn[, 1 ], Tr_wid = Tr_bvn[, 2 ])
# transform back to simplex
Tr_splx <- bvn_to_simplex (
cbind (df$ Tr_loc, df$ Tr_wid))
df$ Tr_loc_splx <- Tr_splx[,1 ] + .5 * Tr_splx[,2 ]
df$ Tr_wid_splx <- Tr_splx[,2 ]
```
### Marginal Distributions
```{r}
df %>%
ggplot () +
geom_histogram (
aes (Tr_loc_splx),
fill = "lightblue" ,
alpha = 1 ,
binwidth = .005
) +
geom_histogram (aes (Tr_wid_splx),
fill = "purple" ,
alpha = .3 ,
binwidth = .005 ) +
scale_x_continuous (
limits = c (0 , 1 ),
labels = c ("0" , ".25" , ".5" , ".75" , "1" ),
expand = expansion ()
) +
scale_y_continuous (limits = c (0 , NA ), expand = expansion ()) +
labs (x = "Marginal Locations / Widths" , y = "Density" ) +
theme_pubr ()
```
```{r}
df %>%
ggplot () +
geom_histogram (
aes (Tr_loc),
fill = "lightblue" ,
alpha = 1 ,
binwidth = .005
) +
geom_histogram (aes (Tr_wid),
fill = "purple" ,
alpha = .3 ,
binwidth = .005 , ) +
scale_x_continuous (
limits = c (NA , NA ),
expand = expansion ()
) +
scale_y_continuous (limits = c (0 , NA ), expand = expansion ()) +
labs (x = "Marginal Locations / Widths" , y = "Density" ) +
theme_pubr ()
```
### Joint Distribution
### Ternary
```{r}
df %>%
ggplot () +
geom_abline (intercept = 0 , slope = 2 ) +
geom_abline (intercept = 2 , slope = - 2 ) +
geom_point (
aes (x = Tr_loc_splx, y = Tr_wid_splx),
size = 1 ,
alpha = .1 ,
shape = 16
)+
scale_x_continuous (
limits = c (0 , 1 ),
labels = c ("0" , ".25" , ".5" , ".75" , "1" ),
expand = expansion ()
) +
scale_y_continuous (
limits = c (0 , 1 ),
labels = c ("0" , ".25" , ".5" , ".75" , "1" ),
expand = expansion ()
) +
theme_pubr ()
```
#### Unbounded Scale
```{r}
df %>%
ggplot () +
geom_point (
aes (x = Tr_loc, y = Tr_wid),
size = 1 ,
alpha = .1 ,
shape = 16
) +
theme_pubr ()
```
## Lambda (Item Discernibility)
### sigma_lambda
```{r}
df$ sigma_lambda <- exp (rnorm (1e5 , log (.5 ), .5 ))
#median(lambda$sigma_lambda)
df %>%
ggplot () +
geom_histogram (aes (x = sigma_lambda),
fill = "lightblue" ,
binwidth = .01 ) +
scale_x_continuous (limits = c (0 , 2.5 ), expand = expansion ()) +
scale_y_continuous (limits = c (0 , NA ), expand = expansion ()) +
theme_pubr ()
```
### lambda
```{r}
df$ lambda_loc <- exp (rnorm (1e5 , 0 , df$ sigma_lambda))
df$ lambda_wid <- exp (rnorm (1e5 , 0 , df$ sigma_lambda))
df %>%
ggplot () +
geom_histogram (aes (x = lambda_loc),
fill = "lightblue" ,
binwidth = .01 ) +
scale_x_continuous (limits = c (0 , 5 ), expand = expansion ()) +
scale_y_continuous (limits = c (0 , NA ), expand = expansion ()) +
theme_pubr ()
```
## E (Person Proficiency)
### mu_E
Then narrow Student t prior ensures that the population proficiency is regularized to 1.
```{r}
df$ mu_E <- rnorm (1e5 )
df %>%
ggplot () +
geom_histogram (aes (x = exp (mu_E)),
fill = "lightblue" ,
binwidth = .01 ) +
scale_x_continuous (limits = c (0 , 5 ), expand = expansion ()) +
scale_y_continuous (limits = c (0 , NA ), expand = expansion ()) +
theme_pubr ()
```
### sigma_E
```{r}
df$ sigma_E <- exp (rnorm (1e5 , log (.5 ), .5 ))
df %>%
ggplot () +
geom_histogram (aes (x = sigma_E),
fill = "lightblue" ,
binwidth = .01 ) +
scale_x_continuous (limits = c (0 , 2.5 ), expand = expansion ()) +
scale_y_continuous (limits = c (0 , NA ), expand = expansion ()) +
theme_pubr ()
```
### E
```{r}
df$ E_loc <- exp (rnorm (1e5 , mean = df$ mu_E, sd = df$ sigma_E))
df$ E_wid <- exp (rnorm (1e5 , mean = df$ mu_E, sd = df$ sigma_E))
df %>%
ggplot () +
geom_histogram (aes (x = E_loc),
fill = "lightblue" ,
binwidth = .01 ) +
scale_x_continuous (limits = c (0 , 5 ), expand = expansion ()) +
scale_y_continuous (limits = c (0 , NA ), expand = expansion ()) +
theme_pubr ()
```
### Lambda and E joint
```{r}
df %>%
ggplot () +
geom_histogram (aes (x = 1 / lambda_loc / E_loc),
fill = "lightblue" ,
binwidth = .01 ) +
scale_x_continuous (limits = c (0 , 5 ), expand = expansion ()) +
scale_y_continuous (limits = c (0 , NA ), expand = expansion ()) +
labs (x = "SD of Latent Appraisal" , y = "Density" ) +
theme_pubr ()
```
## Omega (Residual Correlation)
```{r}
beta_param <- 2
df$ omega <- rbeta (1e5 , shape1 = beta_param, shape2 = beta_param)* 2-1
df %>%
ggplot () +
geom_histogram (aes (x = omega),
fill = "lightblue" ,
binwidth = .01 ) +
scale_x_continuous (limits = c (- 1 , 1 ), expand = expansion ()) +
scale_y_continuous (limits = c (0 , NA ), expand = expansion ()) +
theme_pubr ()
```
## Latent Appraisals
### Marginal Appraisals
```{r}
A <- rbvnorm (1e5 ,
mean1 = df$ Tr_loc,
mean2 = df$ Tr_wid,
sd1 = 1 / df$ lambda_loc / df$ E_loc,
sd2 = 1 / df$ lambda_wid / df$ E_wid,
cor = df$ omega)
df$ A_loc <- rnorm (1e5 , mean = df$ Tr_loc, sd = 1 / df$ lambda_loc / df$ E_loc)
df$ A_wid <- rnorm (1e5 , mean = df$ Tr_wid, sd = 1 / df$ lambda_wid / df$ E_wid)
A_splx <- bvn_to_simplex (df[,c ("A_loc" , "A_wid" )])
df$ A_loc_splx <- A_splx[,1 ] + .5 * A_splx[,2 ]
df$ A_wid_splx <- A_splx[,2 ]
df %>%
ggplot () +
geom_density (aes (A_loc_splx),
fill = "lightblue" ,
alpha = .3 ,
shape = 16 ) +
geom_density (aes (A_wid_splx),
fill = "purple" ,
alpha = .3 ,
shape = 16 ) +
scale_x_continuous (limits = c (0 , 1 ), expand = expansion ()) +
scale_y_continuous (limits = c (0 , NA ), expand = expansion ()) +
labs (x = "Appraisal Locations / Widths" , y = "Density" ) +
theme_pubr ()
```
## Joint Appraisals
```{r}
df %>%
ggplot () +
geom_abline (intercept = 0 , slope = 2 ) +
geom_abline (intercept = 2 , slope = - 2 ) +
geom_point (
aes (x = A_loc_splx, y = A_wid_splx),
size = 1 ,
alpha = .1 ,
shape = 16
)+
scale_x_continuous (
limits = c (0 , 1 ),
labels = c ("0" , ".25" , ".5" , ".75" , "1" ),
expand = expansion ()
) +
scale_y_continuous (
limits = c (0 , 1 ),
labels = c ("0" , ".25" , ".5" , ".75" , "1" ),
expand = expansion ()
) +
theme_pubr ()
```
# Biases
## Shifting Bias
### sigma_b
```{r}
df$ sigma_b <- exp (rnorm (1e5 , log (.5 ), 1 ))
df %>%
ggplot () +
geom_histogram (aes (x = sigma_b),
fill = "lightblue" ,
binwidth = .01 ) +
scale_x_continuous (limits = c (0 , 5 ), expand = expansion ()) +
scale_y_continuous (limits = c (0 , NA ), expand = expansion ()) +
theme_pubr ()
```
### b
```{r}
df$ b_loc <- rnorm (1e5 , mean = 0 , sd = df$ sigma_b)
df$ b_wid <- rnorm (1e5 , mean = 0 , sd = df$ sigma_b)
df %>%
ggplot () +
geom_histogram (aes (x = b_loc),
fill = "lightblue" ,
binwidth = .01 ) +
scale_x_continuous (limits = c (- 5 , 5 ), expand = expansion ()) +
scale_y_continuous (limits = c (0 , NA ), expand = expansion ()) +
theme_pubr ()
```
```{r}
b_splx <- bvn_to_simplex (df[,c ("b_loc" , "b_wid" )])
df$ b_loc_splx <- b_splx[,1 ] + .5 * b_splx[,2 ]
df$ b_wid_splx <- b_splx[,2 ]
df %>%
ggplot () +
geom_abline (intercept = 0 , slope = 2 ) +
geom_abline (intercept = 2 , slope = - 2 ) +
geom_point (
aes (x = b_loc_splx, y = b_wid_splx),
size = 1 ,
alpha = .1 ,
shape = 16
)+
scale_x_continuous (
limits = c (0 , 1 ),
labels = c ("0" , ".25" , ".5" , ".75" , "1" ),
expand = expansion ()
) +
scale_y_continuous (
limits = c (0 , 1 ),
labels = c ("0" , ".25" , ".5" , ".75" , "1" ),
expand = expansion ()
) +
theme_pubr ()
```
## a (Scaling Bias)
### sigma_a
```{r}
df$ sigma_a <- exp (rnorm (1e5 , log (.5 ), .5 ))
df %>%
ggplot () +
geom_histogram (aes (x = sigma_a),
fill = "lightblue" ,
binwidth = .01 ) +
scale_x_continuous (limits = c (0 , 2.5 ), expand = expansion ()) +
scale_y_continuous (limits = c (0 , NA ), expand = expansion ()) +
theme_pubr ()
```
### a
```{r}
df$ a_loc <- exp (rnorm (1e5 , mean = 0 , sd = df$ sigma_a))
df %>%
ggplot () +
geom_histogram (aes (x = a_loc),
fill = "lightblue" ,
binwidth = .01 ) +
scale_x_continuous (limits = c (0 , 5 ), expand = expansion ()) +
scale_y_continuous (limits = c (0 , NA ), expand = expansion ()) +
theme_pubr ()
```
# Responses
## Marginal Responses
```{r}
df$ Y_loc <- df$ A_loc * df$ a_loc + df$ b_loc
df$ Y_wid <- df$ A_wid + df$ b_wid
Y_splx <- bvn_to_simplex (df[,c ("Y_loc" , "Y_wid" )])
df$ Y_loc_splx <- Y_splx[,1 ] + .5 * Y_splx[,2 ]
df$ Y_wid_splx <- Y_splx[,2 ]
df %>%
ggplot () +
geom_histogram (
aes (Y_loc_splx),
fill = "lightblue" ,
alpha = 1 ,
binwidth = .005
) +
geom_histogram (aes (Y_wid_splx),
fill = "purple" ,
alpha = .3 ,
binwidth = .005 ) +
scale_x_continuous (
limits = c (0 , 1 ),
labels = c ("0" , ".25" , ".5" , ".75" , "1" ),
expand = expansion ()
) +
scale_y_continuous (limits = c (0 , NA ), expand = expansion ()) +
labs (x = "Marginal Locations / Widths" , y = "Density" ) +
theme_pubr ()
```
## Joint Responses
```{r}
df %>%
ggplot () +
geom_abline (intercept = 0 , slope = 2 ) +
geom_abline (intercept = 2 , slope = - 2 ) +
geom_point (
aes (x = Y_loc_splx, y = Y_wid_splx),
size = 1 ,
alpha = .1 ,
shape = 16
)+
scale_x_continuous (
limits = c (0 , 1 ),
labels = c ("0" , ".25" , ".5" , ".75" , "1" ),
expand = expansion ()
) +
scale_y_continuous (
limits = c (0 , 1 ),
labels = c ("0" , ".25" , ".5" , ".75" , "1" ),
expand = expansion ()
) +
theme_pubr ()
```
***
# Recovery Check
## Define Parameter Bounds for Variances
```{r}
# compute a benchmark for the mean and SD of the parameters
mean_benchmark <- simplex_to_bvn (c (.4 , .2 , .4 ), type = "ilr" )
sd_benchmark_loc <- simplex_to_bvn (c (.98 , .01 , .01 ), type = "ilr" )
sd_benchmark_wid <- simplex_to_bvn (c (.495 , .01 , .495 ), type = "ilr" )
# mean for Tr_loc
mu_Tr_loc <- mean_benchmark[1 ]
# mean for Tr_wid
mu_Tr_wid <- mean_benchmark[2 ]
# SD forTr_loc
sigma_Tr_loc <- sd_benchmark_loc[1 ] / 4
# SD Tr_wid
sigma_Tr_wid <- abs (sd_benchmark_wid[2 ] - mean_benchmark[2 ]) / 4
# SDs for other parameters
sigma_lambda_E_loc <- .3
sigma_lambda_E_wid <- .3
sigma_a_loc <- .3
sigma_b_loc <- sigma_Tr_loc / 3
sigma_b_wid <- sigma_Tr_wid / 3
```
## Generate Data for One Replication
```{r}
n_respondents <- 50
n_items <- 20
sim_data <-
generate_itm_data_sim_study (
n_respondents = n_respondents,
n_items = n_items,
mu_Tr_loc = mu_Tr_loc,
mu_Tr_wid = mu_Tr_wid,
sigma_Tr_loc = sigma_Tr_loc,
sigma_Tr_wid = sigma_Tr_wid,
sigma_lambda_E_loc = sigma_lambda_E_loc,
sigma_lambda_E_wid = sigma_lambda_E_wid,
sigma_a_loc = sigma_a_loc,
sigma_b_loc = sigma_b_loc,
sigma_b_wid = sigma_b_wid,
omega_beta = 3
)
responses <- sim_data$ responses
```
### Plot Data
```{r}
samples <- gather_values (
lower = responses$ x_L,
upper = responses$ x_U,
item_id = responses$ jj,
step_size = 0.01
)
samples %>%
ggplot (aes (x = samples)) +
geom_density (fill = "purple" , alpha = .5 ) +
facet_wrap (~ item_id) +
scale_x_continuous (
limits = c (0 , 1 ),
labels = c ("0" , ".25" , ".5" , ".75" , "1" ),
expand = expansion ()
) +
theme_pubr () +
theme (plot.margin = margin (.1 , .5 , .1 , .1 , "cm" ),
panel.grid.major = element_line ())
```
### Plot True Intervals as Error Bars
```{r}
data.frame (
idx = 1 : n_items,
lower = sim_data$ parameters$ Tr_L,
upper = sim_data$ parameters$ Tr_U
) %>%
ggplot (aes (y = idx)) +
geom_errorbarh (aes (xmin = lower, xmax = upper), height = 0.3 ) +
scale_x_continuous (
limits = c (0 , 1 ),
labels = seq (0 , 1 , .25 ),
expand = expansion ()
) +
labs (x = "Latent Truth" ,
y = "Item" ) +
theme_pubr () +
theme (plot.margin = margin (.1 , .5 , .1 , .1 , "cm" ),
panel.grid.major = element_line ())
```
### Plot True Intervals as Ternary Plot
```{r}
data.frame (
idx = 1 : n_items,
Tr_L = sim_data$ parameters$ Tr_L,
Tr_U = sim_data$ parameters$ Tr_U) %>%
mutate (
Tr_loc = Tr_L + (Tr_U - Tr_L) / 2 ,
Tr_wid = Tr_U - Tr_L
) %>%
ggplot () +
geom_abline (intercept = 0 , slope = 2 ) +
geom_abline (intercept = 2 , slope = - 2 ) +
geom_point (
aes (x = Tr_loc, y = Tr_wid),
size = 4 ,
alpha = .5 ,
shape = 16
)+
scale_x_continuous (
limits = c (0 , 1 ),
labels = c ("0" , ".25" , ".5" , ".75" , "1" ),
expand = expansion ()
) +
scale_y_continuous (
limits = c (0 , 1 ),
labels = c ("0" , ".25" , ".5" , ".75" , "1" ),
expand = expansion ()
) +
theme_pubr ()
```
## Fit Model
### Stan Data
```{r}
### Stan data declaration
I <- length (unique (responses$ ii))
J <- length (unique (responses$ jj))
N <- nrow (responses)
ii <- responses$ ii
jj <- responses$ jj
nn <- c (1 : N)
Y_splx <- cbind (responses$ x_splx_1, responses$ x_splx_2, responses$ x_splx_3) %>%
as.matrix ()
#Y_splx[is.na(Y_splx)]
# sanity check: all rows sum to 1
# check <- apply(Y_splx, 1, sum) %>% as.data.frame()
# table(check)
## stan data list
stan_data <- list (
I = I,
J = J,
N = N,
ii = ii,
jj = jj,
nn = nn,
Y_splx = Y_splx,
link = 1
)
```
### Compile Model
```{r}
# Choose model to fit
model_name <- "itm_simulation_v2_beta"
# Compile model
model <-
cmdstanr:: cmdstan_model (
stan_file = here ("src" , "models" , paste0 (model_name, ".stan" )),
pedantic = TRUE ,
quiet = FALSE
)
```
### Run Sampler
```{r eval=TRUE, message=FALSE, warning=FALSE}
# number of MCMC chains
n_chains <- 4
# Run sampler
fit_prior_check <- model$sample(
data = stan_data,
seed = 2023,
chains = n_chains,
parallel_chains = n_chains,
iter_warmup = 500,
iter_sampling = 500,
refresh = 500,
adapt_delta = .99,
init = .1
)
# save fit_prior_check
fit_prior_check$save_object(
file = here("fits", paste0(model_name, "_fit_prior_check.rds")))
```
```{r}
# load fit
fit_prior_check <- readRDS (
file = here ("fits" , paste0 (model_name, "_fit_prior_check.rds" )))
```
***
### Get Estimates
```{r}
estimates_summary <- fit_prior_check$ summary ()
```
### Sampler Diagnostics
```{r}
fit_prior_check$ diagnostic_summary ()
```
```{r}
# sampler diagnostics
fit_prior_check$ sampler_diagnostics (format = "df" ) %>%
psych:: describe (quant = c (.05 ,.95 ),) %>%
round (2 ) %>%
as.data.frame () %>%
dplyr:: select (median, min, Q0.05 , Q0.95 , max) %>%
.[- c (7 : 9 ),]
```
```{r}
# convergence diagnostics
convergence_summary <-
fit_prior_check$ draws (format = "df" ) %>%
summarise_draws (.x = ., "rhat" , "ess_bulk" , "ess_tail" ) %>%
remove_missing () %>%
select (- variable) %>%
psych:: describe (., quant = c (.05 , .95 )) %>%
as.data.frame () %>%
select (mean, median, Q0.05 , Q0.95 , min, max)
convergence_summary %>% round (3 )
```
### Effective sample size (ESS) & Rhat Plots
```{r message=FALSE, warning=FALSE}
# color scheme
color_scheme_set(scheme = "purple")
# Effective sample sizes
plot_neff <-
mcmc_neff_hist(bayesplot::neff_ratio(fit_prior_check), binwidth = .01) +
labs(title = "A") +
guides(color = FALSE, fill = FALSE) +
theme(
legend.text = element_blank(),
legend.key = element_blank(),
title = element_text(size = 16, face = "bold")
)
# Rhat
plot_rhat <-
bayesplot::mcmc_rhat_hist(bayesplot::rhat(fit_prior_check)) +
labs(title = "B") +
guides(color = FALSE, fill = FALSE) +
theme(
legend.text = element_blank(),
legend.key = element_blank(),
title = element_text(size = 16, face = "bold")
) +
yaxis_text(on = TRUE)
# Combined plot
plot_diagnostics <- gridExtra::grid.arrange(plot_neff, plot_rhat, ncol = 2)
```
## Parameter Plots
```{r}
# color scheme
color_scheme_set (scheme = "purple" )
true_parameters <- sim_data$ parameters
```
### Person Competence: Location
```{r }
plot_E_loc <-
mcmc_recover_scatter(x = fit_prior_check$draws("E_loc"), true = true_parameters$E_loc) +
labs(
subtitle = "Person Competence: Location",
)
plot_E_loc
```
### Person Competence: Width
```{r }
plot_E_wid <-
mcmc_recover_scatter(x = fit_prior_check$draws("E_wid"), true = true_parameters$E_wid) +
labs(subtitle = "Person Competence: Width")
plot_E_wid
```
### Person Scaling Bias: Location
```{r }
plot_a_loc <-
mcmc_recover_scatter(x = fit_prior_check$draws("a_loc"), true = true_parameters$a_loc) +
labs(subtitle = "Person Scaling Bias: Location") +
theme(
text = element_text(
family = "serif",
size = 16,
colour = 1
),
plot.title = element_text(size = 20,
face = "bold"),
plot.margin = unit(c(.1, .1, .1, .1), "cm"),
axis.title = element_text(size = 16, color = 1),
axis.title.x = element_text(margin = margin(t = 5)),
axis.title.y = element_text(margin = margin(r = 5)),
axis.text = element_text(size = 12, colour = "black"),
axis.text.x = element_text(margin = margin(t = 3))
)
plot_a_loc
```
### Person Shifting Bias: Location
```{r }
plot_b_loc <-
mcmc_recover_scatter(x = fit_prior_check$draws("b_loc"), true = true_parameters$b_loc) +
labs(subtitle = "Person Shifting Bias: Location")
plot_b_loc
```
### Person Shifting Bias: Width
```{r }
plot_b_wid <-
mcmc_recover_scatter(x = fit_prior_check$draws("b_wid"), true = true_parameters$b_wid) +
labs(subtitle = "Person Shifting Bias: Width")
plot_b_wid
```
### Item Discernability: Location
```{r }
plot_lambda_loc <-
mcmc_recover_scatter(x = fit_prior_check$draws("lambda_loc"),
true = true_parameters$lambda_loc)
plot_lambda_loc
```
### Item Discernability: Width
```{r }
# plot
plot_lambda_wid <-
mcmc_recover_scatter(x = fit_prior_check$draws("lambda_wid"),
true = true_parameters$lambda_wid) +
labs(subtitle = "Item Discernability: Width")
plot_lambda_wid
```
### Residual Correlation
```{r }
plot_omega <-
mcmc_recover_scatter(x = fit_prior_check$draws("omega"),
true = true_parameters$omega) +
xlim(-1, 1) +
ylim(-1, 1) +
labs(subtitle = "Residual Correlation")
plot_omega
```
### Latent Truth Intervals
```{r }
# plot
plot_Tr_loc <-
mcmc_recover_scatter(x = fit_prior_check$draws("Tr_loc"),
true = true_parameters$Tr_loc) +
labs(subtitle = "Item True Location")
plot_Tr_loc
```
```{r }
# plot
plot_Tr_wid <-
mcmc_recover_scatter(x = fit_prior_check$draws("Tr_wid"),
true = true_parameters$Tr_wid) +
labs(subtitle = "Item True Width")
plot_Tr_wid
```
```{r}
# ITM estimates
latent_truth_est_ilr <- data.frame (
idx = 1 : J,
type = "estimated" ,
Tr_L_est = fit_prior_check$ summary ("Tr_splx" ) %>%
as.data.frame () %>%
pull (median) %>%
.[1 : J],
Tr_wid_est = fit_prior_check$ summary ("Tr_splx" ) %>%
as.data.frame () %>%
pull (median) %>%
.[(J + 1 ): (J * 2 )]
) %>%
mutate (Tr_loc_est = Tr_L_est + Tr_wid_est / 2 ,
Tr_U_est = Tr_L_est + Tr_wid_est)
# compute means in the unbounded space
means_ilr <-
cbind (jj, simplex_to_bvn (Y_splx)) %>%
as.data.frame () %>%
dplyr:: group_by (jj) %>%
dplyr:: summarise (across (everything (),mean, na.rm = TRUE )) %>%
dplyr:: ungroup () %>%
dplyr:: select (- jj) %>%
bvn_to_simplex () %>%
mutate (
idx = 1 : J,
type = "simple_mean" ,
Tr_L_sm = x_1,
Tr_U_sm = 1 - x_3
) %>%
select (- c (x_1, x_2, x_3))
# true parameters
latent_truth_true <- data.frame (
idx = 1 : J,
type = "true" ,
Tr_L_true = true_parameters$ Tr_L,
Tr_U_true = true_parameters$ Tr_U
)
df_interval_plot_ilr <- full_join (latent_truth_est_ilr, latent_truth_true) %>%
full_join (means_ilr) %>%
mutate (type = factor (type, levels = c ("estimated" , "true" )))
```
#### Interval Plot
```{r message=FALSE, warning=FALSE}
cols <- c("True" = "grey70","ITM (ILR)" = "red", "Mean" = "blue")
df_interval_plot_ilr %>%
ggplot(aes(y = idx)) +
geom_errorbarh(
aes(xmin = Tr_L_true, xmax = Tr_U_true, col = "True"),
height = 0,
linewidth = 5
) +
geom_errorbarh(
aes(xmin = Tr_L_sm, xmax = Tr_U_sm, col = "Mean"),
height = 0,
linewidth = 3,
alpha = .4
) +
geom_errorbarh(
aes(xmin = Tr_L_est, xmax = Tr_U_est, col = "ITM (ILR)"),
height = 0,
linewidth = 2,
alpha = .5
) +
scale_x_continuous(
limits = c(0, 1),
labels = seq(0, 1, .25),
expand = expansion()
) +
scale_y_discrete(
expand = expansion(add = 1)
) +
scale_color_manual(values = cols) +
labs(x = "Latent Truth", y = "Item") +
theme_pubr() +
theme(plot.margin = margin(.1, .5, .1, .1, "cm"),
panel.grid.major = element_line())
```
# Prior vs. Posterior
## True Intervals
```{r}
latent_truth <- data.frame (
Tr_loc_splx_post = fit_prior_check$ draws ("Tr_loc_splx" ) %>%
unlist () %>%
as.vector (),
Tr_wid_splx_post = fit_prior_check$ draws ("Tr_wid_splx" ) %>%
unlist () %>%
as.vector ()
) %>%
mutate (
jj = factor (rep (1 : J, each = nrow (.) / J)),
Tr_wid_splx_prior = rbeta (nrow (.), 1.2 , 3 ),
Tr_loc_splx_prior = (1 - Tr_wid_splx_prior) * rbeta (nrow (.), 1 , 1 ) + Tr_wid_splx_prior / 2
)
Tr_bvn <- simplex_to_bvn (
cbind (
latent_truth$ Tr_loc_splx_prior - .5 * latent_truth$ Tr_wid_splx_prior,
latent_truth$ Tr_wid_splx_prior,
1 - (
latent_truth$ Tr_loc_splx_prior + .5 * latent_truth$ Tr_wid_splx_prior
)
)
)
latent_truth <- cbind (latent_truth,
Tr_loc_prior = Tr_bvn[, 1 ], Tr_wid_prior = Tr_bvn[, 2 ])
latent_truth %>%
ggplot () +
geom_abline (intercept = 0 ,
slope = 2 ,
alpha = .7 ) +
geom_abline (intercept = 2 ,
slope = - 2 ,
alpha = .7 ) +
geom_density_2d_filled (
aes (x = Tr_loc_splx_prior,
y = Tr_wid_splx_prior),
binwidth = .05 ,
alpha = .3 , ) +
geom_point (
aes (x = Tr_loc_splx_post, y = Tr_wid_splx_post, col = jj),
size = .3 ,
alpha = .3 ,
shape = 16
) +
scale_x_continuous (
limits = c (0 , 1 ),
labels = c ("0" , ".25" , ".5" , ".75" , "1" ),
expand = expansion ()
) +
scale_y_continuous (
limits = c (0 , 1 ),
labels = c ("0" , ".25" , ".5" , ".75" , "1" ),
expand = expansion ()
) +
labs (x = "Location" , y = "Width" ) +
# supress legends
guides (col = FALSE , fill = FALSE ) +
theme_pubr ()
```
## Hyper-Priors
### sigma_lambda: Location
```{r}
colors <- c ("prior" = "lightblue" , "posterior" = "purple" )
sigma_lambda_loc <- data.frame (posterior = fit_prior_check$ draws ("sigma_J[1]" , format = "list" ) %>% unlist () %>% as.vector ())
sigma_lambda_loc$ prior <- exp (rnorm (nrow (sigma_lambda_loc), log (.5 ), .5 ))
sigma_lambda_loc %>%
ggplot () +
geom_density (aes (prior), fill = "lightblue" , alpha = 1 ) +
geom_density (aes (exp (posterior * .5 + log (.5 ))), fill = "purple" , alpha = .3 ) +
scale_x_continuous (limits = c (NA , 3 ), expand = expansion ()) +
scale_y_continuous (limits = c (0 , NA ), expand = expansion ()) +
labs (x = "Prior / Posterior" , y = "Density" ) +
scale_color_manual (values = colors) +
theme_pubr ()
```
### sigma_lambda: Width
```{r}
sigma_lambda_wid <- data.frame (
posterior = fit_prior_check$ draws ("sigma_J[2]" , format = "list" ) %>% unlist () %>% as.vector ()
)
sigma_lambda_wid$ prior <- exp (rnorm (nrow (sigma_lambda_loc), log (.5 ), .5 ))
sigma_lambda_wid %>%
ggplot () +
geom_density (
aes (prior),
fill = "lightblue" ,
alpha = 1 ) +
geom_density (aes (exp (posterior* .5 + log (.5 ))),
fill = "purple" ,
alpha = .3
) +
scale_x_continuous (
limits = c (NA , 3 ),
expand = expansion ()
) +
scale_y_continuous (limits = c (0 , NA ), expand = expansion ()) +
labs (x = "Prior / Posterior" , y = "Density" ) +
scale_color_manual (values = colors) +
theme_pubr ()
```
### sigma_E: Location
```{r}
sigma_E_loc <- data.frame (
posterior = fit_prior_check$ draws ("sigma_I[1]" , format = "list" ) %>% unlist () %>% as.vector ()
)
sigma_E_loc$ prior <- exp (rnorm (nrow (sigma_lambda_loc), log (.5 ), .5 ))
sigma_E_loc %>%
ggplot () +
geom_density (
aes (prior),
fill = "lightblue" ,
alpha = 1 ) +
geom_density (aes (exp (posterior* .5 + log (.5 ))),
fill = "purple" ,
alpha = .3
) +
scale_x_continuous (
limits = c (NA , 3 ),
expand = expansion ()
) +
scale_y_continuous (limits = c (0 , NA ), expand = expansion ()) +
labs (x = "Prior / Posterior" , y = "Density" ) +
scale_color_manual (values = colors) +
theme_pubr ()
```
### sigma_E: Width
```{r}
sigma_E_wid <- data.frame (
posterior = fit_prior_check$ draws ("sigma_I[2]" , format = "list" ) %>% unlist () %>% as.vector ()
)
sigma_E_wid$ prior <- exp (rnorm (nrow (sigma_lambda_loc), log (.5 ), .5 ))
sigma_E_wid %>%
ggplot () +
geom_density (
aes (prior),
fill = "lightblue" ,
alpha = 1 ) +
geom_density (aes (exp (posterior* .5 + log (.5 ))),
fill = "purple" ,
alpha = .3
) +
scale_x_continuous (
limits = c (NA , 3 ),
expand = expansion ()
) +
scale_y_continuous (limits = c (0 , NA ), expand = expansion ()) +
labs (x = "Prior / Posterior" , y = "Density" ) +
scale_color_manual (values = colors) +
theme_pubr ()
```
### mu_E: Location
```{r}
mu_E_loc <- data.frame (
posterior = fit_prior_check$ draws ("mu_I[1]" , format = "list" ) %>% unlist () %>% as.vector ()
)
mu_E_loc$ prior <- rnorm (nrow (sigma_lambda_loc))
mu_E_loc %>%
ggplot () +
geom_density (
aes (prior),
fill = "lightblue" ,
alpha = 1 ) +
geom_density (aes (posterior),
fill = "purple" ,
alpha = .3
) +
scale_x_continuous (
limits = c (- 4 , 4 ),
expand = expansion ()
) +
scale_y_continuous (limits = c (0 , NA ), expand = expansion ()) +
labs (x = "Prior / Posterior" , y = "Density" ) +
scale_color_manual (values = colors) +
theme_pubr ()
```
### mu_E: Width
```{r}
mu_E_wid <- data.frame (
posterior = fit_prior_check$ draws ("mu_I[2]" , format = "list" ) %>% unlist () %>% as.vector ()
)
mu_E_wid$ prior <- rnorm (nrow (sigma_lambda_loc))
mu_E_wid %>%
ggplot () +
geom_density (
aes (prior),
fill = "lightblue" ,
alpha = 1 ) +
geom_density (aes (posterior),
fill = "purple" ,
alpha = .3
) +
scale_x_continuous (
limits = c (- 4 , 4 ),
expand = expansion ()
) +
scale_y_continuous (limits = c (0 , NA ), expand = expansion ()) +
labs (x = "Prior / Posterior" , y = "Density" ) +
scale_color_manual (values = colors) +
theme_pubr ()
```